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Abstract 

The third order optical response of bosonic excitons is calculated using the Green's function solution of the Nonlinear Exciton 
Equations (NEE) which establish a quasiparticle-scattering mechanism for optical nonlinearities. Both time ordered and non 
ordered forms of the response function which represent time and frequency domain techniques, respectively, are derived. New 
components of the response tensor are predicted for isotropic ensembles of periodic chiral structures to first order in the optical 
wavevector. The nonlocal nonlinear response function is calculated in momentum space, where the finite exciton-exciton 
interaction length greatly reduces the computational effort. Applications are made to coupled anharmonic vibrations in the 
amide I infrared band of peptides. Chirally-sensitive and non sensitive signals for a helices and antiparallel f3 sheets are 
compared. 
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I. INTRODUCTION 



Calculating the nonlinear optical response of large molecules using conventional sum-over-states expressions is a 
challenging task since the number of states accessible by multiple quantum excitations increases rapidly with system 
size Q For example, the Frenkel exciton model for M. coupled three level systems has M. one-exciton states 
and M x (Ai + l)/2 two-exciton states. Diagonalizingthe two-exciton hamiltonian is the bottleneck in numerical 



simulations. The Nonlinear Exciton Equations (NEE) |3j, U, 10 offer an alternative exact method for computing the 
exact nonlinear optical response of systems whose Hamiltonian conserve the number of excitons. The NEE establish 
an exciton scattering mechanism and provide a practical algorithm for computing the third order optical response, 
totally avoiding the calculation of two-exciton states. The NEE were first developed by Spano and Mukamel |a,flQ 
and applied for four-wave mixing of coupled two-level and three-level III, ISl 111 , 121 molecules. The local- 

field-approximation was generalized by adding two-exciton variables which properly account for double quantum 
resonances. Additional dynamical variables have subsequently been added, allowing to describe population transport 
via the Redfield equation for the reduced exciton density matrix [lil Q| . The NEE were further extended to particles 
with arbitrary commutation relations, and to Wannier excitons in semiconductors P. ll5ill^ |. We have recently studied 
the frequency-domain third order susceptibility of isotropic exciton systems and calculated the leading terms of its 
tensor elements going beyond the dipole approximation [17j . In this paper we apply these results to the time domain 
response function in a collinear field configuration, which leads to the strongest signals. 

The NEE provide a unified treatment of vibrational and electronic excitons. The Frenkel excitons corresponding 
to electronic excitations are Paulions whose non-Boson operator statistics results in nonlinear hard-core re pul sive 
interactions between excitons. The NEE were originally derived using a Bosonization procedure |fjl7tl8l ll3Ulflllfll2(l ] 
whereby the Pauli exciton operators are replaced by Boson operators and a repulsive potential is added to the 
hamiltonian. This results in a soft-core boson model with a finite anharmonicity A; the original hard-core boson 
hamiltonian is recovered by setting A — > oo. 

This article focuses on vibrational excitons which are intrinsically bosons QE3Q and the finite anharmonicities 
serve as the source of nonlinearities. Third order impulsive optical techniques performed with linearly polarized light 
(figure^) are commonly used to probe ultrafast processes in isotropic systems. The corresponding nonlinear response 
function S^, aVaVl is a fourth rank tensor ,21] which relates the third order polarization, P^(r,i), at position r and 
time t to the optical field E(r',t'); Vi = x,y,z are the polarization components of the i'th laser pulse in the lab 
frame. In general the response function has 3^ tensor elements. For isotropic systems in the dipole approximation, 
only those elements with 5 ViV;i 5 V2Vl , 8 ViU2 8 V3Vl and 5 ViVl 5 V3V2 (e.g. xxxx, and zzyy) survive rotational averaging, and 
three linearly-independent elements xxyy, xyxy and xyyx, are necessary to describe the response for an arbitrary 
pulse polarization configuration j^l^- However, other elements with an odd number of repeating indices (such as 
xxxy, zzxy), which vanish in this approximation due to isotropic symmetry, appear when the dipole approximation 
is relaxed. We have recently shown that these elements are chir ally-sensitive i e appear only in chiral molecules 
("handed" systems which are distinct from their mirror images |23l l24j ~). and vanish in racemates, equal mixtures of 
molecules with opposite sense of chirality. 
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The dipole approximation implies that the optical electric field is uniform across the molecule, thus, its phase factor 
exp(zkr), where k is the wavevector, does not affect the response and can be set to zero. However, the wavevector does 
play an important role in the spectroscopy of chiral molecules, where the variation of the phase induces new tensor 
components of the response function to first order in the wavevector, which vanish in the dipole approximation. The 
various tensor elements can be probed directly using different time-domain techniques which control the sequences of 
optical interactions. 

The difference in absorption of left- and right-handed circularly polarized light [3, 13, l^fil known as circular 
Dichroism (CD) is the simplest example of a wavevector-induced signal and is related to the s2y elements of the 
linear response tensor (when the field propagates along z). The technique is sensitive to molecular structure and is 
extensively used for protein structure determination both in the UV (180 - 220 nm n — it* and it — tt* transitions) 
and the IR (1000 - 3500 cm"'' which covers most of the amide vibrational bands) j^, I30I I31I l^ . Pattern- 
recognition and decomposition algorithms are used to distinguish between a-helical and /3-sheet formations using 
electronic QQQ and vibrational CD Q. 

Polypeptides often have almost translationally-invariant secondary structures (helices, sheets, strands) which form 
periodic arrays of localized vibrations. Analyzing secondary structure motifs, thus, provides specific information which 
can be used for studying globular proteins with different secondary structures. Periodicity can be used to greatly 
reduce the computational effort, as is the case for electronic excitations in molecular crystals and semiconductor 
superlattices ^fj, [3- [ill I41I . Due to different translational properties of one-exciton and two-exciton 
states in the Frenkel exciton model, multi-exciton resonances cannot be generally calculated analytically even for 
infinite periodic systems. The NEE only require the one-exciton states, and yield closed expressions for infinite 
periodic systems, where translational symmetry helps reduce the problem size even further. Relaxing the dipole 
approximation using the NEE, is straightforward. The CD spectra of molecular aggregates modelled as collections of 
coupled electric dipolcs were calculated to first order in wavevectors |43j. This model has been applied to biological 

nn 

light harvesting antennae and cylindrical aggregates JTJj, |45J . By extending this procedure to the nonlinear response 
we obtain the complete set of tensor elements for the response function of infinite periodic systems. 

In section[n]we present the time-domain expressions for the third order optical response. The hamiltonian and the 
NEE for vibrational excitons are given in section ITTT1 and the third order response function is derived in section llVl 
Two techniques for probing one-exciton and two-exciton resonances are discussed in section The two-dimensional 
infrared frequency correlation signals of two typical structural motifs of polypeptides: the one-dimensional a helix and 
the two-dimensional antiparallel (3 sheet (figure^) in the amide I region (1550 - 1750 cm~^) are compared. The results 
are discussed in section IVII The linear absorption is calculated in appendix Rotational averages for isotropic 
systems are given in appendix [BJ The nonlinear signals are calculated in appendices El (four- wave mixing in various 
phase matching directions) and [D] (pump-probe). The frequency domain signal is given in app endix El and the exciton 
scattering matrix for an infinite periodic system is presented in appendix IfI 
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II. TIME-DOMAIN OPTICAL RESPONSE OF EXCITONS 



The optical response of molecules is determined by the induced polarization, which serves as the source in the 

(n) 

Maxwell equations for the generated signal field. The response functions S{,J +1 v n ...vi are system property-tensors 
which allow to calculate the induced polarization for an arbitrary incoming pulse configuration |2| : 

pw +i ( x „ +1 )= f dx «-/ dXi 

Here E is the optical electric field vector and v — x, y, z denote its Cartesian components in the lab frame; x = (r, t) 
is the space-time vector and J dx = J dr J dt, where the r integration runs over the molecular volume, while the t 
integration is from — oo to +00. x represent the times and coordinates of the interaction with the optical pulses. 

Using causality S^ n ' vanishes unless t\ . . . t n precede t n+ \ for all n. We focus on time-domain experiments where 
the interaction sequence is controlled by short optical pulses and the time-ordered response function is finite for 
t n +i > t n > ... > t2 > t\ and vanish otherwise as shown in figure Et- Alternatively we can require that the response 
function to be symmetric with respect to permutation olvjitj, j — 1,2,3. This non-time-ordered response function, 
which is useful for frequency domain techniques is given in appendix [E] 

The linear polarization is given by: 

p ll ] (x2) = E / dxiSgU^.xiJlMx!), (2) 

where S$ Vl (xa, xi) is the nonlocal linear response function. is responsible for linear absorption (local response; 
see appendix |A"|) and CD (non local response). Four wave mixing (4WM) signals are described by the third order 
polarization 

Pl 3 1 ) (x 4 ) = I dx 3 / dx 2 f d Xl 

J J J 

S i%^ ( x 4, x 3) x 2 , *i)E„ 3 (x 3 )E U3 (x 2 )E ul ( Xl ). (3) 

We assume three short well-separated incoming pulses: 

1 3 

E v (x) = -^2El s) (t-t s )exp{ik s r-ioj s (t-t s )+i<f> s ) + c.c, (4) 

s=l 

where Ei s \t — t s ) is the (real) envelope of pulse s centered at t s , with wavevector k s , carrier frequency lu s and phase <p s . 
When the pulses are much shorter than the relevant molecular timescale, E v (t — t s ) in eq. (2) can be approximated as 
e[ , S(t — t s ). Since the pulses are longer than the optical periods we must then invoke the rotating wave ap prox imation 



46J. In this 



- KWA and only retain S^ 3 ' terms resonant with optical fields which fall within the pulse bandwidth oj s 
approximation the field amplitude in the frequency domain is taken to have a rectangular shape centered at w s with 
width uj s . We next assume that the pulses are temporally well separated and ordered, i. e. the first pulse 
interacts at t\, followed by E^ at t<i and the last pulse is E^ at t 3 . The third order polarization is then given by 

P^)(x4) = exp(±^ 3 ± i<f> 2 ± i<t>i) 
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x ^3 E S^ BBtfl (k4t4 > ±k3*3,±k2*2,±kit 1 )^^£;W, (5) 

V3U2V1 

where x = (k, i) ; tj and kj now coincide with the central pulse times and wavevectors and we have applied the 
space-time Fourier transform F(k,u>) = J dt J drexp(iu>t + ilcr)F(r,t). 

When the system is initially at equilibrium, the response function is time translationally invariant and only depends 
on the positive time intervals T s = t s+ i — t s (see figure |2J). Space translational invariance for isotropic systems implies 
k 4 = =pk 3 Tk 2 Tki. We thus denote P^(-k 4 ,< 4 ) = P^(T 3 ,T 2 ,T 1 ) with k s = -k 4 . There are four independent 
signal wavevectors kg: kf = — ki + k 2 + k3, k// = +ki — k 2 + k3, km = +k 4 + k 2 k3 and k/y = +ki + k 2 + k3. 



III. THE NEE FOR VIBRATIONAL EXCITONS 



The amide vibrations of polypeptides can be modeled as N coupled anharmonic vibrational modes localized at the 
peptide bonds and described by the exciton hamiltonian: 

m^n 

H ^ ^mB m B m -\- ^ Jm^nB\ n B n 
m. ran 

B\ n B\B m ,B n , - J rfrP(r) ■ E(r, r). (6) 



mn,ra n 



The creation, B^, and annihilation, B m , operators for mode m satisfy the boson [B m ,Bl] = S mn commutation 
relations. The first two terms represent the free-boson harmonic hamiltonian where e m is the frequency of mode to, 
and the quadratic intermode coupling, J m ,m is calculated in the Heitler-London approximation where off resonant 
B m W n and B m B n terms are neglected. The third term represents a quartic anharmonicity. We assume a pairwise 
anharmonic interaction, U mn ,m'n' = A 4 '" {5 m m'5 n n' + S mn >S nm i), where A mm is the on-site anharmonicity of the 
overtone band and A nm is the intermode anharmonicity of the combination band. These anharmonicities constitute 
the exciton-exciton scattering potential. For A = the hamiltonian describes free bosons, which is a linear system 
whose nonlinear response vanishes identically [if . The fourth term in the hamiltonian represents the interaction with 
the optical field E(r, t), where 

P(r) = ^5(r-r ro ) Mm ( J Bt 1 + B ro ) (7) 

m 

is the polarization operator and \x m is the transition dipolc moment of mode m located at r m ; a vector with components 

(MmiAi Mm)- 

The expectation value of the polarization operator which describes the vibrational response to the optical field will 
be calculated by solving the NEE 0, 0, Q] ■ This hierarchy of equations of motion for exciton variables may be exactly 
truncated order by order in the field since the molecular hamiltonian conserves the number of excitons, and the optical 
field creates or annihilates one exciton at a time. When pure-dephasing is neglected, the only required variables for 
the third order response are B m = (B m ) (one-cxciton) and Y mn = (B m B n ) (two-exciton) and the NEE read QjI^]: 

— i I" ^ ] h m , n B n — fi m {j) ~ ^ ] Vml'm'n' Bi,Y m i n i 7 (8) 

n I'm'n' 
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~ i fa + X! Kn'mV + V mn ,m'n') Y m'n' = fim(r)B n + p, n (T)B m . (9) 

m'n' 

here /i m ,„ = <5 m ,„£ m + J m , n (l - <S m ,„) is an effective one-exciton hamiltonian, ft^2 , m '„' = &m',mh n ,n' + 5 n>n -h m ,m' is 
a free-two-exciton hamiltonian, y mnjm - n ' = f7 m „, m v + U nm , m > n > = -f^^mm'^n' + 5mn'<W') is the anharmonicity 
matrix and /} m (r) = /x m • E{v m ,r) (for the actual times we use r instead of t which were defined as time-ordered). 
The polarization is given by the expectation value of eq. J7J): 

P(r, r) = S (* ~ r m ) Mm B m (r) + c.c. (10) 

m 

The nonlinearities in these equations originate from the anharmonicity: as indicated earlier, for V = eq. JSJ is 
linear, and the nonlinear response vanishes (the two-exciton variable can be exactly factorized as Y rnn = B m B n in eq. 

© )■ 

The evolution of a single-exciton created by an impulsive excitation, fill' (t), is described by the one-exciton Green's 
function G(r): 

Bi! ) W=EG m , m( (r)B< 1 , ) (0), (11) 

m' 

which satisfies the equation: 

dG m Ar) +i J2h m , n ,G n , !n (T)=6(r). (12) 
This equation can be solved using the one-exciton eigenenergies, E^, and eigenvectors, ip^ m - 

^ h m>n l/>en = E^Vfm- (13) 

n 

We then get 
where 

k( T ) = 0{T)exp(-iE S T -7 f r), (15) 

7^ is a dephasing rate of the £ exciton state, and the Heavyside step function (9(t) — for r < and 0(r) = 1 for 
r > 0) guarantees causality. 

Similarly we define the two-exciton Green's function Q Y : 

Y mn {r) = ]T Q* ntm , n , (r)Y rn , n , (0), (16) 

m'n' 

whose equation of motion: 

^ ^ ' Q l mn,m"n" Vmn,m"n")G m "n" .m'n' ~ ^( T )- (1^) 
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Q Y is connected to the zero-order noninteracting (V = 0) two-exciton Green's function Q by the Bethe-Salpeter 
equation 



G Y {t) = Q{t)+ f dr' f dnG(r - t')T(t' -t x )Q{t x ), 
Jo Jo 



(18) 



which defines the two exciton scattering matrix T. Both the two-exciton Green's function and the scattering matrix are 
tetradic matrices; like G(t), the scattering matrix is causal as well. Q can be factorized into a product of one-exciton 
Green's functions, Q mn ,m'n'{^) = G m , m '(r)G n) „/(r). 

The frequency domain scattering matrix T(lu) obtained by the Fourier transform, T(uj) = J dtexp(iujt)T(t), is 
calculated in appendix IF1 P. Iri| : 

r(w) = -iV(l + ig(u)V)- 1 . (19) 

where 

6mn,m'n'(w) = £ V^m^'n^f' MV>| m ' V>|'„- > ( 20 ) 

with 

%'M = \ ■ ( 21 ) 

Calculating the scattering matrix requires the inversion of the matrix Z) = 1 + iQ(cj)V whose matrix elements: 

Dmn.iji,^) — &miO~nj ~l~ ^ ^ ^ Qmn.m'n' (^ J )^ / m'n / ,ij ■ ( 22 ) 
m'n' 

The required numerical effort can be reduced considerably for periodic systems and when the short range nature of 
anharmonicities is taken into account |3j. 



IV. THE NONLINEAR OPTICAL RESPONSE: GREEN'S FUNCTION SOLUTION OF THE NEE 



The nonlinear response is calculated by an order-by-order expansion of the NEE variables in the field using the 
exciton Green's functions, where the optical field and the lower order variables serve as the sources. Setting B m = 
B$ + B$ + B$ + ... and Y mn = Y^i + + Y ( H + ... we get to third order Ubj: 



714 V 



2i / dr" / dr' / dr 3 / dr 2 / dr x 

J — OO J — OG J — OO J — OO J — OO 

J2 £ e{T 2 -T X )T nMn{ {T"-T') 

xG n4t< (r 4 - r")Gl, „ 3 (r" - r 3 )G n>2 (r' - r 2 )G„ ; , ril (r' - n) 

X An3( T 3)/ i ™ 2 (T-2)/i™ 1 (Tl), (23) 
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The r' and t" variables denote the times of the first and the last exciton-exciton interaction respectively, as shown 
in figure Et- The third order polarization is finally obtained from eqs. IjlOjl and (|23J) : 

/n /• poo />00 />oo 

dr 3 / dr 2 / dr x I dr 3 / dr 2 / dn 
J J J oc J oc J — CO 

E (MnfnXn, (r4,r 3) r 2 , ri )) 



n4n 3 n 2 ni 



/OO /"CO 
dr" / dr 1 £ ^2-n)r nin ^ nWl (r"-r') 
-OO J — OO III! 



xG„ 4 ,„, (r 4 - r")Gl, „ 3 (r" - r 3 )G„^„ 2 (r' - r 2 )G„- ni (r' - n) 

x^3(r 3 r 3 )£; i/2 (r2T2)£; i , 1 (riTi) +c.c, (24) 



where 



M^™ 1 (r 4 ,r3,r 2 ,r 1 ) = <5(r 4 - r„ 4 )<5(r 3 - r„ 3 )<5(r 2 - r„ 2 )5(r x - r„J 

x^^:, (25) 

"c.c" is the complex conjugate and (...) denotes rotational averaging (appendix [BJ). ft is important to note that 
unlike t±, t 2 and i 3 , the integration variables n, T2 and t 3 do not have any particular time ordering. 

Since the Green's functions are retarded (i. e. they vanish for negative time arguments) only three sequences of 
interaction times contribute to eq. I|24|) : i) t 2 > T\ > t 3 , ii) T2 > T3 > ti, iii) r 3 > T2 > ri. In each case we switch 
to a different set of time-ordered variables: i) £4 = T4, £ 3 = t 2 , t 2 = t\, t\ = t 3 , ii) £4 — T4, i 3 = t 2i t 2 — t 3 , ii = n, 
iii) ^4 = T4, t 3 = r 3 , £2 = r 2 , £i = tj (figure |5Jd). Eq. © together with eq. 124|) then result in the following three 
contributions to the response function: 

^4 ) .. l , 1 ( x 4,--,xi) = 5^ ) .^ 1 (x 4 ,...,xi) 
+S£>„ 1 (x 4) ...,x 1 ) 

+^ 7 . / ) 1 (x 4 ,...,x 1 ) + c.c. (26) 

Each term now corresponds to one particular interaction sequence ( is obtained from (i), S^ 11 ^ - from (ii) and 
S^ 111 ^ - from (iii) ). These are given by: 



Sgl Vl (x 4 , Xl ) = 2i E (M^n, fa, r 3 , r 2 , r x )> 

TI4. . .711 

x/ dr" dr' s V r n / n /„^ n / (r s " - t' s ) 
Jo Jo 1 , , 1 

xG n4<( T s) G n^ 3 ( i « - T' s ')G n ' 2 n 2 (U 2 - r ")G J t l ,„ 1 ^41 ~ t' s ), (27) 

Sg 1 ^ (x 4 , Xl ) - 2» E (M^T^ ( r 4, r 8) r a , r x )> 

7T.4. . .711 

x / dr'J / d<^ E r «S^«;( T "-^) 
Jo Jo i , i i 



71 4 ?lo Tl^Tl^ 

2»2 



xG n4nl (ri)G„, n3 (i43 ~ r>')Gl (t« - t^)G„ w (t« - O. ( 2 «) 
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Sg'X (x 4 , Xl ) = 2* Y, < M S« ( r 4, r a> r 2 , n)) 



X 



£4 2 



E 

n'. ni.ni.ni 



dr'J dr' s Y Fn'Mn'M'-ri). 
o Jo 



xG„ 4nl (^)G£, n3 (i 4 3 - ri)G n>2 (t 42 - OG„i»x(*« - (29) 

Here t' s denotes the time interval between the polarization detection and first exciton scattering event, while t'J denotes 
the interval between the detection and the last exciton scattering event, as shown in figure [5^,. Using these time- 
ordered expressions we can switch to new variables representing the time intervals between interactions tjj = ij — tj 
with i > j: £43 = T3, £42 = T3 + T 2 , ^41 = T3 + T2 + Ti. Thus, all terms in the response function only depend on 
these three positive time intervals. So far, the three terms in the response function merely represent different time 
orderings in the integrations, however, we will shortly see that they represent distinct optical signals. 

The time intervals in eqs. (|27|) ~ H29|I and their relations to the actual interaction times are depicted in figureEl The 
response can be interpreted using figure Et: Let us consider S^ 111 ': three interactions with the optical fields at times 
t\ < t-2 < £3 generate three quasi-particles (two with positive oscillation frequency and one - with negative) whose 
evolution is described by the one-exciton Green's functions G n ' n , , G n ' n „ and , . The two positive frequency 
quasiparticles at n\ and n' 2 are scattered by r n / n ^ n / to n' 3 and n' A . One of the two excitons, n' 3 , corresponds to the 
exciton generated by the third field, n' 3 *— n 3 . The other exciton, n' 4 — » 71,4, generates the optical response. and 
S^ 11 ^, differ by the scattering sequence of events as shown in figure EJ; and may be interpreted similarly. 

Using eq. and invoking the same approximations used in section^ we can relate the wavevectors and times in 
the response functions (eqs. H27JI - 129JI ) to the wavevectors and times of the optical pulses. Assuming that the carrier 
frequencies of all pulses are resonant with the one-exciton manifold, we select the resonances in the response functions 
using the RWA and express the nonlinear polarization for the signal wavevectors k/, k// and k///. To find out the 
wavevector dependence we transform eqs. l|2*?jl - (l2l31) to momentum space using F(k) = J drF(r) exp(ikr). Then: 

PlliT^Tx) = ^ exp(+^ 3 + ifc - #1) 

(T 3 ,T 2 ,Ti)£g> (30) 



l/ 3 V2 V\ 



P k 4 "(T 3 , T 2 , T x ) = i exp(+z0 3 - i<h + #1) 



E ^Uz^T^EifEifEW, (31) 



x 



p£" (T 3 , T 2> Ti) = ^ exp(-# 3 + ifc + %<h) 



x s k "/ 



(T^^T^ElfElfElW (32) 



/_ ^ 1/4 1/3 U2 V\ \ & ' ^ ' ± / V3 V2 U\ 



where we have denoted 



-V4 U3 U2 V 1 



(T 3 , T 2> Ti) = ((ki - k 2 - k 3 )t 4 , k 3 t 3 , k 2 i 2 , -Mi), (33) 
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Vi u 3 u 2 v i 



fill! 



(Ta.Ta.Ti) = ((k 2 - k x - k 3 )i 4 ,k 3 i 3 , -M^kjA), (34) 



(T 3 , T 2 ,Ti) = ((kg - k 2 - ki)t 4> -Ma, k 2 t 2 , Mi). (35) 



W^e, i/iws, /ind i/iai £/ie response function generates the kj signal, S^- 11 ' generates hjj and S^ 111 * 1 generates k///. 
The fourth possible signal P kiv vanishes for the present model since it has no transition dipole connecting the three- 
exciton states with the ground state. The three terms in eqs. I|27 fl -I|29 [) were obtained by a simple bookkeeping of 
time variables. The RWA has connected these terms with the impulsive signals in eqs. (|30|) - (|35|) . Each of the three 
signals is thus given by a single term. 

Equations (|27|) - (|29[) are used in appendix IO to calculate the signals in the eigenstate basis. The frequency domain 
scattering matrix (see appendix [FJ allows to simplify response functions in eqs. 

lEUl-dHl- The three signals are 
given by eqs. IjClljl . (JC14II and (|C17|) . The sequential pump-probe spectrum is calculated in appendix IDl 

By transforming all time variables to the frequency domain, §(f2 3 , f2 2 , f2i) = 
J °° dT 3 r o °° dT 2 J °° f2Ti§(T 3 , T 2 ,T 1 ) exp(^ 3 T 3 + ifl 2 T 2 + ifiiTi), (see appendixEJ equations JSl"© give: 

s^^(ns,n a ,ni) = a £ K:( k * - k ? k 3 )dg*(-k 3 )dg* (-k 2 )d£(-k x )> 

x J &&( n a)^ 1 (-n l )^(n3)r^ 1 £, Ca (n a + e 6 + l ^ 1 )i^ 2 (n 3 + e Si (36) 

fe^^fi,,^) = 2t Kl^ k 3 - k 1 )dg*(-k 3 )dg(-k 2 )d^*(-k 1 )) 

«4-..«l 

X^x( n 2)^x(0l)^4(^3)r^ 2 , 6 e 1 (n 3 + E £2 +i 7£2 )I £3£i (0 3 + E S2 +i 7fo ), (37) 

(ns,^, no = 2* y <c( k 3- k 2- k i)rfg(- k 3)4r(- k 2) d sr(- k o> 

«4-..«l 

x % (ni)/ f4 (fi 3 )/|, (fia - o 3 )[r^ 3>?26 (o 2 )% 6 (o 2 ) 

- r «4« 3 ,«2?i (O3 + e ?3 + i7« 3 )^2«i (O3 + E 4a + i7£ s )], (38) 
where we have defined the exciton transition dipole for state £ 

^(k)-^e ikr - M ^^ m (39) 

JTi 

and 

'«<") = 5T^W5' (40) 

n 

The application of the Green's functions expressions for periodic infinite systems is straightforward |17j : the sum- 
mations over one-exciton eigenstates are replaced by summations over different Davydov exciton bands at zero mo- 
mentum. The scattering matrix of infinite systems which involves all possible momenta of different Davydov bands 
is given in appendix IfI 
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V. APPLICATION TO THE AMIDE I BAND OF PEPTIDES 



We have applied the present theory to the amide I vibrations of two ideal structural motifs of polypeptides: a 
helices (one dimensional) with 18 amide residues in the unit cell and antiparallel (3 sheets (two dimensional) with 4 
residues per unit cell. The structural and coupling parameters were reported earlier [itIL^ . The anharmonicity A m „ 
is local (A = -16 cm - ' for m = n and zero otherwise) and the same dephasing rate, 75 = 3 cm - ' was assumed for all 
excitons. We used 100 cells in each dimension to calculate the scattering matrix with periodic boundary conditions 
(eq. (US )■ 

The linear absorption of both motifs presented in figure shows two absorption peaks for the a helix and the (3 
sheet. For the a helix the longitudinal transition - along the helix axis gives a peak at 1642 cm - ', and transverse - in 
the plane perpendicular to the helical axis gives a weaker peak at 1661 cm"' [43 ■ For the f3 sheet the main peak at 
1632 cm - ' is horizontal (parallel to the sheet surface). The other weak peak at 1700 cm - ' consists of two transitions 
[4^: vertical (perpendicular to the sheet surface transition) at 1707 cm"' and horizontal at 1699 cm - '. 

Signals were calculated for a collinear field configuration where, ki, k2 and k3, propagate along z. We further 
assume all fields to have the same carrier frequency so that |kj| = uj s /c for all kj (j — 1, 2, 3) where c is the speed of 
light. The absolute magnitudes of the S** V3V2Ul (Cls, T 2 — 0, fii) and ^tl" 3V2Vl (^3, ^2, ?i = 0) signals were computed 
using eqs. I|C11|) and i|C17() . These are one-sided Fourier transforms of the time domain k/ and km signals. 

The &£l V3 v2 Vl signals for the a helix are shown in the left column in figure The xxxx component is finite in the 
dipole approximation and shows one major peak associated with the longitudinal transition [47j and a much weaker 
20 cm - ' blue shifted transverse transition. The crosspeaks between these two peaks are symmetric with respect to the 
diagonal and have roughly the same amplitudes. This signal resembles our previous calculation of large (90 residue) 
helices Q]> which could be treated as infinite helices with small edge defects. 

This xxxx pattern is changed in the chirally-sensitive xxxy component. The diagonal peaks are suppressed. This 
can be explained by noting that the distance between sites enters the signal amplitude: diagonal peaks have no such 
distance factor since they originate from interactions with the same mode. The crosspeaks are induced by interactions 
of different modes and depend on their distance. Thus, in the chiral signal the crosspeaks are amplified compared to 
the diagonal peaks. Additionally, the crosspeak pattern is asymmetric since one of the crosspeaks (below the diagonal) 
dominates. The other, above the diagonal, crosspeak is much weaker. The remaining two chiral components xxyx 
and xyxx also suppress the diagonal peaks, while show two crosspeaks 

The antiparallel (3 sheet signals displayed in the right column in figure 01 show a similar pattern, however only 
one horizontal diagonal peak (as shown in ^^1) 1S visible in xxxx with very weak crosspeaks shifted to much higher 
frequencies compared to the diagonal. Since there are four sites per unit cell, we generally expect four optical 
transitions (Davydov components). This could result in four diagonal peaks and six crosspeaks in both sides of the 
diagonal if all transitions are correlated. However, xxxx is dominated by a single strong horizontal transition. Similar 
to the helix, this peak structure changes in the chiral component, xxxy, which shows much stronger crosspeaks than 
the diagonal peaks. The strongest crosspeak is well separated from the strongest xxxx peak. Like the helix, the xxyx 
and xyxx configurations are both similar and show crosspeaks above and below the diagonal. 
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We next turn to the hjjj technique. ^'" U2l/ shows two-exciton resonances along L02 (double quantum coherence); 
the resonances along W3 originate both from one- to two-exciton resonances (e to / in figure ^) and from one-exciton 
to ground state resonances (0 to e). The xxxx component for the a helix presented in figure [5] shows one major 
peak and a weaker peak originating from the same two-exciton resonance. The main peak of the chiral components 
xxxy, xxyx, xyxx is shifted compared to xxxx, indicating that different two-exciton resonances make the strongest 
contribution to the signal. The xxxy and xxyx signals are very similar and show two strong peaks unlike xyxx which 
only shows one major peak. The xxxx component of the /3 sheet displayed in figure has a similar structure to the 
a helix. The chiral xxxy component shows two peaks at higher two-exciton resonance comparing to xxxx. xxyx, is 
very similar, xyxx shows one major peak at 103 — 1700 cm" . 



VI. DISCUSSION 



In the electric dipole approximation, the linear absorption of isotropic systems is related to the diagonal tensor 
elements of the linear susceptibility |l7l |43j . Nonlinear signals calculated within this approximation then provide 
a limited window onto the optical responses of isotropic ensembles through three independent chirally-non-sensitive 
tensor components of the third order response function: xxyy, xyyx, xyxy. The signal propagation direction is 
determined by phase-matching. Going beyond the dipole approximation, we found three additional components for a 
collinear laser configuration: xxxy, xxyx, xyxx. Noncollincar configurations (which can also satisfy phase matching) 
lead to six nonzero elements (the three additional components are: (z)xyzz, (z)xzzy and (z)xzyz, where the first 
index (z) denotes the wavevector, and the other four are polarization components). These are chirally-sensitive and 
show a different pattern in the correlation plots, which reflects the polarization properties of optical transitions. For 
example, from figures 01 and [S] we see that the strongest peaks originate from correlations between perpendicular 
transitions: longitudinal - transverse in the a helix and horizontal - vertical in the (3 sheet (see figure 

The differences between the chiral components can be deduced by following the interaction and evolution sequences: 
in S k/ , fix corresponds to the free evolution after the first interaction. Subsequently two simultaneous interactions 
take place. Therefore xxyx and xyxx reflect a similar excitation pattern: there is one interaction first, followed by 
two interactions with the perpendicular polarization, xxxy is qualitatively different: the first interaction is followed 
by two interactions with parallel polarizations. This is why, xxxy and xxyx show a different peak pattern. For S k/ " 
there are two simultaneous interactions first, followed by one interaction. Thus, xxxy and xxyx give very similar 
signals since in both cases the first two interactions are with perpendicular fields. The xyxx component is qualitatively 
different since the first two interactions are with parallel fields. 

The signal further carries information regarding the redistribution of excitonic amplitudes between states with 
different polarization properties. For instance, ^y XX originates from exciton amplitude transfer between a state with 
x polarization to a state with y polarization during T2. S^'J y is qualitatively different since the evolution during 
T2 corresponds to x polarized exciton amplitude transfer to another x exciton. S^ xxy again originates from exciton 
amplitude transfer during T\ between states with perpendicular (y — > x) polarizations; xyxx does not require such 
exciton evolution. 
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In the present work we only included homogeneous line broadening caused by fast fre quen cy fluctuations. Slow 
fluctuations result in inhomogeneous broadening which should affect the 2D lineshapes |43, |4±1- The ideal peak 
patterns predicted for periodic structures can be used for structure decomposition and parameter determination 
of real structures. Dynamical information on exciton evolution can be obtained by varying T2 in § k/ (^3, T2, fii). 
Population transport may be incorporated using the theory developed in [lil I2 1| , but this goes beyond the scope of 
this paper. 

We next compare the exciton scattering mechanism offered by the NEE with the more conventional picture of 
transitions among eigenstates, as described by double-sided Feynman diagrams. The relation between the two for kj, 
k// and kj// is shown in figure |SJ For k/ the exciton coherence during T\ in both pictures is represented by the same 
Green's function. However, during T3 the Feynman diagrams show two independent coherence evolutions, only one of 
them involves the two-exciton states. The scattering representation is more complex: during T3 the evolution consists 
of scattering + two free evolutions. This is natural in the molecular basis set where the exciton pathways can be 
followed in real space: the short range exciton scattering is then followed by free evolution as the excitons separate. 
The role of the distance parameter is much less obvious in the eigenstate picture, kj/ is similar, except the density 
matrix evolution during the first and the last intervals has the same frequency sign. In k/77 two excitons are created 
by the initial excitation. Thus the scattering can occur at any time. In the sum over states representation this is given 
as independent evolutions of two-exciton coherences during T± . Both representations of the response are equivalent 
as long as pure dephasing is neglected. The NEE and the corresponding response become more complicated by pure 
dephasing which induces incoherent population transport and requires additional dynamic variables I13I I2 lj | . 

The response function in eq. @ is not required to be time ordered, thus, we can define it to be symmetric with 
respect to permutation of different time arguments. This choice is useful for overlapping optical pulses and frequency 
domain response, where time ordering of interactions is not enforced. Frequency domain signals can be expressed 
using the third order susceptibility given in terms of the exciton scattering matrix as shown in appendix [E] The 
susceptibility is then directly related to non time ordered response function. Expressions for k/, k// and k/// can be 
derived using the RWA. 

Our theory may be applied to other periodic (infinite) as well as to non periodic (finite) systems. Periodicity and 
cyclic boundary conditions result in exciton band structure where only zero momenta exciton states are active. This 
applies as long as the optical wavevector is small compared with the exciton momentum. Exciton-exciton interactions 
are the source of anharmonicities and are described u sing quasi- particle scattering. This description may also be 




applied to semiconductors and quantum superlattices |16t l40t 1 12L |5fJ. When the exciton coherence size becomes 
comparable to the optical wavelength the theory needs to be extended to include polariton effects . 
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Appendix A: LINEAR ABSORPTION OF EXCITONS 



In this appendix we connect the linear absorption with the response function. The linear response function and the 
linear polarization were defined by eq. J5J. The linear absorption can be calculated using the linear susceptibility %W 
|it|. The linear response, especially for vibrations, is often measured in the time domain using ultrashort laser pulses: 
Fourier transform of the time dependent induced polarization (free induction decay) gives the Fourier Transform 
Infrared (FTIR) spectrum. 

We assume a single excitation pulse with envelope E„ (t — to) centered at to with the carrier frequency luq and 
wavevector ko: 

E v (x) = ^8 VUo E^(t- 1 ) exp(ik r - iu (t - to)) + c.c, (Al) 

The linear response function of the excitonic system is given in terms of the one-exciton Green's functions |ri|: 

^W i (x 2 ,x 1 ) =i^(d^(k 2 )4 1 *(-k 1 ))7 ? (t 2 -tO +c.c.' (A2) 
( 

n 

Within the dipole approximation rotational averaging (eq. 1^2(1 *) gives (d^ 2 (k2)<i^ 1 *(— ki)) w i^ 2iyi |<ij| 2 [13. Sub- 
stituting eqs. IjAlj! and l|A2J) into eq. J5J and using assumptions for the field given in section[n]we can integrate J cfoci 
which leads to: 

P„„(k 2 t 2 ) = l -Ei°J J2 \di\ 2 [S(k 2 + ko)7 e (t 2 - fi) - S(k 2 - k )^(t 2 - h)}, (A3) 
S 

where <5(k) accounts for the translational invariance of an isotropic system. This expression shows that _P„ (k 2 t 2 ) 
only depends on the delay T\ = t 2 — ti and can, thus, be denoted i- > t/() (k 2 t 2 ) = P Vo (k 2 , T\). The first term in eq. 
(|A3|) in real space is proportional to e lk ° r and, thus, describes forward propagation of the induced polarization. 

The absorption spectrum is defined as the imaginary part of Fourier transform of the linear polarization with respect 
to T\. <ja{u) — Ira ^ dT\ exp(iwTi)i : ' 1/0 (k 2 , Ti), while the real part of this integral describes dispersive lineshapes. 
By eliminating the prefactor {1/Q)E^} '(5(k 2 +ko) and keeping only terms resonant to the positive frequency we obtain: 

The CD spectrum must be calculated by going beyond the dipole approximation 
Appendix B: THE RESPONSE OF ISOTROPIC ENSEMBLES 

Optical fields, wavevectors, and space coordinates are defined in the lab frame, while the transition dipoles and 
their position vectors are given in the molecular frame. The rotational averaging in eqs. (|C3|) and 1|C5(I . (...), needs to 
be performed over the relative orientation of the two frames in order to calculate the response functions for isotropic 
(randomly oriented) ensembles of molecules [22I : 

(d^(k 4 )dg(k 3 )4f(-k 2 )d^*(-k 1 )) 
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x(ei k 4 r„ 4+ ik 3rn3+ik2r „ 2+ ik 1 r„ 1<> ,3^ /[i ,x )) (B1) 

In the phase- matching directions the microscopic response functions only depend on the relative positions of molecules 
and are independent on the origin of the coordinate system. Since the coordinates r m vary only within one molecule, 
for molecules smaller than the wavelength of light we have kr m <C 1, and the exponential function in the transition 
dipole of eq. (|39|l can be expanded to first order. This leads to 

/ ik 4 r„ 4 +ik3r„,+ik 2 r„ 2 +ikir„ 1 vt v 3 v 2 ,.vi \ 
\ c r J, n i r J, n 3 r J 'n2r"'n 1 ' 

+iJ2KKXn>n>n>Z)- ( B2 ) 
In the exciton basis, eqs. I )B1J ( and l|B2|) give 

(^(k 4 )dg(k 3 )^f(-k 2 )^*(-k 1 )) 

K 
ft 

^kf^dgdgdg), (B3) 

K 

where we have used the fact that ip£,m are real and have defined the transition dipole vector for zero momentum 
exciton state 

d£ = d£(k = 0)=^/C^,m, (B4) 

m 

and the tensor 

dT" = y>m/< i ^, m - (B5) 



Eq. (|B3I) requires fourth and fifth rank rotational averagings. The first term in this equation corresponds to 
the dipole approximation. The remaining terms which contain the wavevector and a coordinate, represent a first 
order correction to that approximation. These terms do not depend on the origin of the molecular frame provided 
k 4 + k 3 + k 2 + ki =0 (i. e. the signal wavevector kg = — kj), which is the phase matching condition. 
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Rotational averaging can be performed using the transformation between the lab and molecular frames [2 



(a^...a^)^((e^.a s )...(e^. ai ))= £ tW,,,,...^^ ...a? 1 , (B6) 

a s ...Qi 

where Tv}...v lt a,...ai = {lv s a a ---lviai) is the average of transformation tensor where l VOL is the cosine of the angle 
between laboratory frame axis v = x,y, z and molecular frame axis a — x, y, z. The necessary averages of ranks four 
and five transformation tensors, which are universal quantities independent of system geometry are given in table ?? 

Q 

Using table ??, we obtain for the rotational averages of the transition dipoles: 

(d^d^d^d" 1 ) - T (4) d ai d a3 d a2 d ai CR7) 

(d^d" 3 d" 2 d u /) = V T<$ „ „ „ „ „ „ „ „ <£ 5,a *d?»d?»c£\ (B8) 

N S4 £3 S2 Si' / / Kl'4"3"2"l,a5«4f3(l2ai £4 {3 {2 SI V ' 

a 5 ...ai 

The remaining averages can be simply obtained by permutation of indices. 



Appendix C: TIME-DOMAIN FWM SIGNALS 

The response functions (eqs. Ij27(l - (|29|) ) can be readily calculated using the single exciton basis (eq. I|13(l ) where 
the number of terms in the response function is considerably reduced. We define: 

/" T 2 r T s 

^453?24"l( T 3,T'2,Tl) = / dT'J ^s r ?4«3:?2?l( r " - T s) 

J Jo 

x Iu (Org (rs - r' 3 )I C2 (r 2 - r'J)^ (n - r'J) (CI) 

with the eigenstate basis scattering matrix: 

^Ui3^2il( T ) = 5Z V , { 4 m 4 V'|3m3rm4)n3ni2"il(' r ) 1 /'i;2ni2V';;i)ni- (C2) 

7714 . . .7711 

By transforming the coordinates to momentum space, we obtain from eqs. 1|27|) - (|29|l : 
S$.. V1 (x 4 ,...,x 4 ) 

= 2i £ K:(k4)dg*(-k3)dg*(-k 2 )dg(k 1 )>J^ eifofo (t4i,*43,*«), (C3) 

*4...£l 

^"Li (x 4 ,...,x 4 ) 

= 2i J] (d^k^d^C-ks^k,^ (C4) 
61-6 



1/4... r/i V 



X 4 , ...,X 4 ) 



2i £ (dg(k 4 )dg(k 3 )dg*(-k 2 )d^*(-k 1 )) l 7 e4£3S2ei (t 43 ,t 42 ,t 41 ). (C5) 

«4-4l 
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The wavevector-dependence enters through the transition dipoles. These expressions need to be rotationally averaged 
for isotropic ensembles, as described in appendix 151 
Equations fC^ - fCE^ 

can be conveniently expressed in terms of the frequency domain scattering matrix (eqs. I|19|) 
and (|F1I) ). These expressions then involve triple integrals (two with respect to time and one with respect to the 
scattering matrix frequency). However, the integration limits in eq. (|C1|) are controlled by multiple 9(t) functions 
coming from the Green's functions. Taking these into account and using the frequency domain scattering matrix (eq. 
()19Jl ). the integrals over the exponential functions can be considerably simplified. The integration limits (—00, r 2 ) for 
t'J and (0,t 3 ) for t' s always hold for S 1 , S 11 and S 111 in eqs. (jO^-CSl- However, for S 1 (eq. (O ) and S 11 (eq. 
(|C4|) ) (but not for eq. (|G5() ) other integration limits (0,T2) for t' s can also be used. These considerations allow to 
calculate time integrals: 

sW. 1/1 (x 4 ,...,x 1 ) = 2 J2 <da(k*)dg*(-k3)dg*(-k 2 )d£(ki)> 

?4-..«l 



xi| i(i21 )/ £(i32 )%(t3 2 ) J ^ ^ _ % + f(<y »_ 7fa) (06) 

^ 1 (x 4 ,...,x 1 ) = 2 J] (dg(k 4 )dg*(-k 3 )dg(k 2 )d^*(-k 1 )) 

?4-..«l 

/" 4. f*43) - e- luJti3 I* (t 43 ) 

Xl 6 (t 21 )7| 2 (t 32 )%(t32) / gr g4g2foa (a,)% a (a,) -i^ — ^ (C7) 

^") 1 (x 4 ,...,x 1 ) = 2 J] (^(k 4 )^(k 3 )^*(-k 2 )^*(-k 1 )) 

xkM / ^r fa& 6,fa H , ' , ■ . y ■ (08) 

7 27r W - E «4 - E ?3 + H7f 4 ~ 7? 3 ) 

Equations l|C6l) - (|C8|) constitute our most general expressions for the time domain nonlocal third order responses in 
momentum space. The wavevectors and times in these expressions correspond to the interaction events with the fields. 
All three response functions depend only on time delays betwee different interactions. 

We next turn to on the three signals defined by eqs. (|30|) - (|32J) . The polarization evolution is commonly transformed 
to the frequency-domain where the signal spectra are displayed. The technique k/ is known as photon echo. According 
to the Feynman diagrams (figure EJ) the system is transferred to a coherence after the first interaction. The second 
interaction leaves the system either in a population or in a coherence between two excitonic states. We hold the second 
delay time, T 2 (often referred to as population or waiting time) fixed. Then the population and coherence evolution 
can be probed. The third interaction creates coherences either between the ground and one-exciton states or between 
one- and two-exciton states. To display the signal we perform a double Fourier transform with respect to the first 
and third time delays: T\ — ► SX\ and T3 — » SI3 in eqs. (|30[1 -H32 [1 . The signal also depends on the directions of optical 
wavevectors. Thus, the Fourier transform of eq. Il.'idll can be performed analytically since the Green's functions given 
by eq. (|15|l are simple exponential functions. We finally obtain 

P%(Sl 3 , T 2 , Sl!) = 1 + - vh) E S »Uw (fis.Ta, Sl 1 )E^E^E^. (09) 
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In terms of the scattering matrix and the one-exciton Green's functions, eq. (|C6|I gives: 

/OO J 
arr e4fofe fo(a;)T fofo HJ|' 1 (a; - n 3 ). (CIO) 
-oo z7r 

The frequency integration can be calculated analytically as follows. We note that Y(uj)I{lo) ~ £(cj — 2e + 2i7) _1 is a 
two-exciton Green's function with poles in negative immaginary half-plane (bars indicate averages), while 1^ (cj — i7 3 ) = 
— i(w — ^3 — Egj — ^^j) -1 has one pole in positive immaginary half-plane: lu p — fi 3 + + 17^. In this case we can 
use the Cauchy integral formula by adding integration contour at the positive imaginary half-plane. Then 

S&^ns.Ta.nO = 2z £ (d£(k! -k 2 - k 3 )dg* (-ka)dg* (-k 2 )d£(-k x )) 

x/ s * 1 (T 2 )/ e2 (T 2 )7| 1 (-n 1 )/ e4 (n 3 )r e4Sl ^ 2 (03 + E 5l +i7 Sl )% e2 (n 3 + E Sl +< 7ei ). (Cll) 

Similar to k/, k// (which does not show an echo) can be defined with the same Fourier transform and the same 
time delays: 

^'(ns.Ta.fii) = lexp(+^ 3 -#2 +^1) ]T S^^(n 3 ,T 2) 00^^^. (C12) 

WS VI V\ 

In terms of the scattering matrix and the Green's functions we get from eq. l|G7ll : 

S^' l/3V2Ul (n 3 ,T 2 ,Q 1 ) = 2i ]T (^(ka-ka-kO^^-kaJdgC-kaJd^^-ki)) 

?4-..«l 

/CO ^ 
^U^M^sMlU^ - ^) ( C13 ) 
-00 Z7r 

and after integration over frequency 

S&^tfia.Ta.nx) = 2z ]T <d£(k a - k 3 - k x )dg* (-k 3 )dg(-k 2 )d£*(-kO) 

«4...£l 

4(T 2 )/ 4l (r 2 )/ ?1 (fi 1 )/ 44 (r! 3 )r^ 2 ^ 1 (r! 3 + e C2 +i 7& )^ 3?1 (n 3 + e C2 + i 7& ). (cu) 



x 



This expression is very similar to k/. 

For k///, two interactions with the delay T± create two-exciton coherence with ki + k 2 . By performing the Fourier 
transforms with respect to the second and third time delays: T 2 — > 2 and T 3 — > f2 3 we can observe the two-exciton 
coherence along i7 2 and mixed, 0-to-one and one-to-two, coherences along f2 3 : 

^"(Jls.na.Ti) = ^expHfc + + ]T S^"^^ (0 3 , ^2, T X )E^E^E^, (C15) 

which is obtained from eq. l|C8jl as 

tlllu^^Ti) =2zJ2 (C( k 3 - k 2 - ki)dg(-k 3 )dg*(-k 2 )dg*(-k 1 )) 



■ k , 1 T, ) f,-M >, > j ^Tc 4iM (u)Il (cu - n 3 )6(n 2 - to). (C16) 
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The function 6(10) = i(u> + 27') 1 is taken in the limit 7^ 3> 7' > 0. Analytic integration over frequency involves 
more terms. Again, T(uj)T(lu) ~ i(uj — 2e + 2i^f)~ 1 is a two-exciton Green's function with poles in the negative 
imaginary half-plane, while 1^ (lu — Q, 3 ) = — i(uj — — E^ — i7f x ) has one pole in positive immaginary half-plane: 
Lu p = ^3 + E^ + ij^ . 9(fl 2 — to) = — i(u — 0,2 — 27') _1 has another pole with positive immaginary part: lu P2 — fl 2 +iy. 
In this case Cauchy integration results into two terms and 

§%li„ 2U1 (n 3 ,n 2 ,T 1 ) = 2i K:(k 3 -k 2 -k 1 )^(-k 3 )^f (-k 2 )dg*(-k x )> 

«4..4l 

xi^Ti)/^)/^ - n 3 )[Te dC „ €!lCl (n 2 )j e!lCl (n a ) 

- r «4«3,?2?i( fi 3 + e 6 + iTfoJ^b^COa + E ?3 + i7 fe )], (C17) 

where we have used 7^ 3> 7'. 



Appendix D: THE SEQUENTIAL PUMP-PROBE SIGNAL 



Sequential pump-probe is an incoherent two-pulse FWM technique commonly used for monitoring excited state 
dynamics. The pump and the probe pulses are assumed well-separated and characterized by their carrier frequencies 
ujp U for pump and uj pr for probe, and their delay r pp . We ignore population transport and assume that the delay be- 
tween pump and probe shorter than the population evolution time. In addition, we assume that the pulse bandwidths 
are much narrower than the exciton bandwidth (which always holds for electronic spectroscopy). 

There are two interactions with the pump pulse (ki = ±k ptl , h 2 = Tkpu , where k pu is the pump pulse wavevector) 
and one with the probe (k3 = ±k pr , where k pr is the probe wavevector). We assume that the pump has various 
polarization components with phases 4> pu in eq. ifljl. which can describe e.g. circular polarization of pump. The 
probe can also have all polarization components with phases <j) pr . The pump-probe experiment measures change of 
absorption of probe with respect to its linear absorption. The absorption of the optical field is given by 

a A = J dx[^^^(x)], (Dl) 

where P„(x) is a third order induced polarization and _E„(x) is a probe field when considering pump-probe experiment. 
Using the response function (eq. @) we obtain for the pump-probe signal: 



<T P p(ci! pr ,T pp ,ui pu ) = — j dt 4 J dt 3 J dt 2 J dti 

E%\U - tprMrKts - t pr )E^\t 2 - t^E^ih - t pu ) 



{ 



9 £1/4 1/3 1/2 Vi ( — kp r ^4, -|-kp r ^3, +kp U ^2 5 ~^pu^l) 



dU 

x ex P [z^v(i 4 - i 3 ) - iu P u(t 2 - h) + + 4%. + <fi* - 4!*)] 

x exp[iuj pr (U - U) + iu> pu (t 2 - h) + + cb P r ~ 4> v p l + <P V P u)\ 

0-0 

dti 
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x exp[-iuj pr (U - t 3 ) ~ iu pu (t 2 - h) + t(+0£ - + <ft* - 0£)] 

x exp[-icj pr (t 4 - t 3 ) + *^ pM (i 2 - h) + i(+<f>£ - 4> v p % - <f>; 2 u + 4&)]} (D2) 

Invoking the RWA and the other assumptions regarding the pulses made in section we find that this signal is a 
combination of the kj and k// techniques: 

x/m{s£ Wl ( Upr , r pp , - Wjw )e*(-*»+*J»+*?.-^ 

(<W PP) Wpu )e i C-«+«-«+«) } (D3) 

where k/ = k pr + k pu — k pM and k// = k pr — k p „ + k p „ . It is possible to measure chiral components of this signal 
using left and right handed circular polarizations so that non chiral components cancel. When the pump bandwidth 
is larger than the exciton bandwidth eq. (|D3|I needs to be integrated over uj pu . 

Appendix E: NON-TIME-ORDERED RESPONSE FUNCTION AND FREQUENCY-DOMAIN FWM SIG- 
NALS 



Using eq. I|24|l we can define a non time ordered response function by: 

^ 4 ... 1/1 (X4,...,X 1 ) = | ]T Y ( M ZZn£Aw^ri)) 

perms n^.-.Tii 

/OO rOO 
dr" dr' Y r„ i „^„ ; (r // -r / ). 
-OO J —OO I I I I 

*m ' wi' <Y\' wi' 



J>\rA 



xG nini (U - r")Gl 3n3 (T" - h)G n>2 {T' - t 2 )G n{ni {T' - h) + c.c, (El) 

where ^2 perm3 denotes three terms in the following permutation: (^3^3, V2*2t2, ^l^iti), (vi^iti, ^3^3, vi^iti) and 
(viYiti, viViti, v 3 Y 3 t 3 ). This form can be also expressed in the eigenstate basis using frequency domain scattering 
matrix. 

Eq. IjEljl can be substituted into eq. © to give a non time ordered representation of the response. This form is 
very convenient for frequency domain four- wave-mixing processes : 

P^{k 4 ,uj 4 ) = Y j ' dk 3 J dLJ 3 J dk 2 J dcu 2 J dki J dwx x 

X^lv 3 v> Vl (- k 4 - ^4; k 3 w 3 , k 2 w 2 , kxu)i)E U3 (k 3 ,uj 3 )E„ 2 (k 2) uJ2)E Vl (ki,o;i), (E2) 



where the susceptibility is obtained from eq. (|E1 



Q: 



xS„3^2^i(~ k4 ' -W4;k 3 ,W3,k 2 ,u;2,ki,u;i) = 27ri(5(w 4 - w 3 - w 2 - wi) 



4 E E K:(k4)dg(-k 3 )dg*(k 2 )d^*(k 1 )) 

perm 3 £ 4 £ 3 £ 2 £l 

xr^ 3 ,6«>2 +wi)/ C4 (w 4 )J^(-w 3 )-f&(^)/c 1 (a;i) +c.c.'. (E3) 

20 



where X) pe ™ 3 denotes a sum over the three permutations: (y 3 k 3 ti; 3 , i/ 2 k 2 w 2 , fikiuJi), (z^ 2 k 2 w 2 , ^3k 3 w 3 , v\\s.iui\) and 
(i/ikiu/i, V2k.<zu)2, VsksUs) (the expression is already symmetric to the permutation of z^ikicji and zy 2 k 2 w 2 ). 

is) 

We consider CW laser fields characterized by their amplitude E), , wavevector k s and optical frequency ui s : 

E„(k, u) = ^-E^e^'SQc + K)S(lo - lj s ) + c.c! (E4) 

We consider the signal in the direction — ki — k2 + k3: 

p( 3 4 )(-k 1 -k 2 + k 3 ^ 1 +c 2 - W3 ) = ^ e-^+^+^ElfEifEW 

V3V2U1 

xXv'lvawifci + k 2 - k 3 , -wi - lo 2 + w 3 ;k 3 , -UJ3, -k 2 ^2, -ki^i), (E5) 



In the RWA this gives 



X^luzwfol + k 2 -k 3 ,-Wi ~u 2 +w 3 ;k 3 ,-w 3 ,-k 2 ,w 2 ,-k 1 ,w 1 ) = 2ni 



xr £4«3,«2«i( ti '2 +^i)/ 44 (wi +w 2 - w 3 )/? 3 (w 3 )i£ 2 (w 2 )%(wi) +c.c.'. (E6) 



This is a three dimensional signal of optical frequencies. Two dimensional sections at L03 = io\ of that signal can be 
used Q- 



Appendix F: THE EXCITON SCATTERING MATRIX FOR PERIODIC STRUCTURES 

In this appendix we calculate the time domain scattering matrix (eq. (jCl|) ) for periodic systems, which can be easily 
extended to non periodic finite systems. Using eq. (|19|l we calculate the scattering matrix in the frequency domain, 
which is directly used in appendix |0 The time domain scattering matrix is then given by the Fourier transform: 

r «4?3?2«l W = j ^ CX P(-*^) r «4«3?2«l( W ) ( F1 ) 

An infinite system is constructed by replicating a unit cell where all vibrational or electronic modes are fixed at 
particular sites inside each cell in some crystal lattice. We will assume a cubic lattice of dimensionality T> with a finite 
number N v cells. To guarantee translational invariance we use cyclic boundary conditions. Each mode is represented 
by a pair of indices, Rm, where R is a position vector of the cell and m is the index of the site within the unit cell; 
the position of the m-th site inside the molecule is given by the vector R + p m where R is the origin of the unit cell 
and p m is the displacement from that origin. We denote the number of sites in the cell by Ai and the lattice constant 
a. Since the system is translationally invariant, the intermode coupling JR m .R'„ = J m! „(R' — R) now depends on 
the distance between cells R' — R and on the site indices of each cell, m and n. When R' = R, J m>n (0) describes 
the coupling of modes inside the cell. J m) „(R' — R) with R' 7^ R defines inter-cell couplings. Similar to this coupling 
the intermode anharmonicity AR mj R'„ = A.,„ ,„(R' — R) now depends on the distance between cells and on the site 
indices of each cell. We note that A miJn (R' — R) with R' 7^ R defines anharmonicity of the combination band, where 
two excitations are localized on different sites. 
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The one-exciton eigenstates £ of a periodic system are characterized by two quantum numbers: the Davydov band 
index A is related to different sites in the unit cell, and the exciton momentum (Bloch wavevector) q. For molecular 
systems much smaller than the optical wavelength only zero momentum, q = 0, exciton states contribute to nonlinear 
optical response. The infinite size only enters into the scattering matrix, where excitons with different momenta can 
be involved in the exciton scattering process. 

The equations for the optical response IjClOjl - (|C16|I can be equivalently used for periodic structures with cyclic 
boundary conditions provided the eigenstates £ are replaced with the periodic system eigenstates A at momentum 
q = 0. The scattering matrix is then considerably simplified: 

rA 4 A 3 A 2 Ai (w) = ^ IpXim^Xsma [ X/ ^r"m 4 m 3 ;r'm 2 m 1 (^>)]'ip\ 2 m 2 'ip\ 1 m 1 , (F2) 
rn l ...m 1 — Z c <r'r"</ e 

where r ri , mi „ ;r2im / jn < (oS) is the mixed space scattering matrix (taken at zero momentum). This scattering matrix is 
given in momentum space with respect to the translational motion of excitons, while it explicitly depends on real 
space coordinates ri and r2, which are the distances between pairs of cells - these are not translationally invariant. 
This scattering matrix is given by jl?! : 

f ri ,m,n;r 2 ,m',n'M = -iA ro „(ri) (^(w))^ m ,n;r a ,m',n' ( F3 ) 

and the matrix 

-^ri ,m,n;r2 ,m' .n' (^) 

— ^ri ,r 2 ^m,m' ^n,n' ~t~ ^Qri ,m.n;r2 ,m' .n' (^-0 ^m'n' (^2 ) : (F4) 

where the mixed space two-exciton Green's function 

6n,tn,n;r 2 ,m',n'(w) = y e~ lq(r2 ~ ri) 5 mi „ ;m ',„< (q, -q, Lj) , (F5) 

q 

involves the sum over two-exciton Green's function of one unit cell with different momenta: 

. < _ _ _ _ 

Sm, n ;m'y(q,-q,w) = 2^ tpx m (q)V'A'n (-q)^AA' (q, -q, uW\ m > (q) V'aw (-<!)• ( F6 ) 

AA' 

We have adopted the following notation. ip\ m is a zero momentum wavefunction of exciton band A obtained from 
J2 m ' ^m,m'(q = 0)$\m> = E^Am, where J m ,m'(q = 0) = J2 r Jm,m'{r) and Ea is the eigenenergy of this state. The 
prime in the sums denotes the summation either over sites within one cell or over different Davydov bands at zero 
momentum, while the sum over r is a sum over cells including r = within the scattering length l c defined by 
A ra „(r > l c ) — 0. All site indices m run within one cell, while is the sum over momenta. For infinite systems this 
sum becomes an integral over all exciton bands. 

The scattering matrix calculation (eq. IjCf j l) is the same for the finite non-periodic system except that r and q 
are set to in eqs. (IF2|) - (IF6|) . site indices m and n then correspond to different modes in the entire system and the 
exciton band indices A are changed into exciton states £. 
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Figure Captions 

Figure 1. (top) Time domain third order experiments: three short laser pulses with wavevectors ki, k2 and k3 generate 
a nonlinear polarization with wavevector kg. (bottom left) Structures of the peptide backbone in a helix and 
antiparallel (3 sheet structures: green - C atoms, blue - N atoms, red - O atoms. C=0 responsible for amide 
I mode (1600-1700 cm~^) are emphasized, (bottom right) Energy level scheme of excitonic system with one 
ground state, a manifold of one-exciton states (e) and a manifold with two-exciton states (/). 

Figure 2. a) Time variables used in calculating the excitonic response. Red peaks indicate laser pulses, blue area cor- 
responds to the exciton-exciton scattering process. t\, t 2 and £3 are the first, second and third interactions 
with laser pulses. t± is the signal generation time. = tj+i — i; are the delay times (always positive) 

between two interactions, r" and t' s are the delay times of the exciton scattering, b) Transformation from the 
non ordered time variables ti, t 2 and T3 to ordered times t\, t 2 and t% in eqs. H24I26J) which defines the three 
different scattering pathways, c) Three scattering pathways: k/ = — ki + k2 + k3 involves scattering of the 
excitons created by k2 and kj, k// = ki — k2 +k3 involves scattering of the excitons created by ki and k3, and 
k/// = ki + k2 — k? involves scattering of the excitons created by ki and k2 . 

Figure 3. Linear absorption of the a helix (top) and the antiparallel /? sheet (bottom). 

Figure 4. Absolute value of &5l V3V2Vl T<x = 0, SI 1 ) signal of the a helix and antiparallel (3 sheet (eq. (jClOl) ). Shown 
are one nonchiral, xxxx, and three chiral, xxxy, xxyx, xyxx, components as indicated. Blue crosses mark the 
crosspeaks of xxxx. 

Figure 5. Absolute value of the S^ll' 3ll Vl (^3, ^2, T\ = 0) signal of the a helix and antiparallel j3 sheet (eq. (jClGIl ). Shown 
are one nonchiral, xxxx, and three chiral, xxxy, xxyx, xyxx, components, as indicated. 

Figure 6. Comparison of the exciton-scattering and the transition-among-eigenstates pictures of the three signals 
§ k '(n 3 ,T 2 ,n 1 ), § k "(ft3,T 2 ,fii) and § k "' (fl 3 , n 2 , Ti). fii, tt 2 and fi 3 are the Fourier transform variables 
conjugate to T\, T 2 and T3. Red dots represent the interaction of the system with the field, blue dots mark the 
exciton scattering space. 
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